home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
basic
/
chaosexe.zip
/
XBASINPO.TRU
< prev
next >
Wrap
Text File
|
1980-01-01
|
5KB
|
145 lines
!PROGRAM XBASINPO(INCARE)
CLEAR
PRINT" ***PENDULUM - BASINS OF ATTRACTION***"
PRINT"This program calculates the average angular velocity for"
PRINT"an array of initial conditions:"
PRINT"velocity (-3,3) and angle (-pi ,pi). If the average angular velocity"
PRINT"from a given initial point is positive then circle is placed at the"
PRINT"point, otherwise the average is negative. This program can also "
PRINT"superpose the appropriate Poincare section on the graph."
PRINT"The program can also save the data sets to two different files."
PRINT
LIBRARY "sglib.trc"
DIM a(1),b(1)
DECLARE DEF accel
INPUT prompt"Input driving force strength: ":g
INPUT prompt"input damping :":q
INPUT Prompt"Input min. and max. time of averaging:":tmin,tmax
INPUT Prompt"Poincare attractor yes (1), no(2) ":Poincare
INPUT prompt"Save data to a file, yes(1), no(2):":sv
IF sv=1 then
PRINT"Basin of attraction File name includes:"
PRINT" 1)First 2 digits for q value"
PRINT" 2)Next needed digits for g value"
PRINT" 3)Last digits as 0's"
INPUT prompt"file name - format ex. 20150000:":filename
INPUT prompt"data file drive a,b,c,etc.?":b$
IF poincare=1 then
PRINT"Superposed Poincare section file name is similar to above except"
PRINT"that the numeral '0' is added, for example, 20015000."
INPUT prompt"Input corresponding Poincare file name:":poinfile
END IF
LET FILENAME$=STR$(FILENAME)
LET POINFILE$=STR$(POINFILE)
END IF
!
CALL PARAMS(W,EPS,TSTEP)
CALL SETXSCALE(-3,3)
CALL SETYSCALE(-3,3)
CALL SETTEXT("PENDULUM - BASINS OF ATTRACTION","INIT. ANGLE","INIT. ANG. VEL.")
CALL RESERVELEGEND
IF sv=1 then
OPEN #1:name b$&":"&filename$,organization record,create newold
ASK #1: FILESIZE length
IF length=0 then SET#1:rECSIZE 10
SET #1: POINTER end
IF poincare=1 then
OPEN #2:name b$&":"&poinfile$,organization record,create newold
ASK #2: FILESIZE length
IF length=0 then SET#2:rECSIZE 10
SET #2:pOINTER end
END IF
END IF
DATA 0,0
CALL DATAGRAPH(A,B,1,0,"WHITE")
CALL gotocanvas
LET phi=0
FOR xint=-pi to pi step .1
FOR vint=-3 to 3 step .15
LET x= xint
LET v=vint
LET t=0
LET s=0
FOR i=1 to 1000000
CALL rk4(x,v,tstep,xnew,vnew,t,w,g,q) ! Take a step of size tstep
LET tshalf=tstep/2
CALL rk4(x,v,tshalf,xnh,vnh,t,w,g,q) !Take two half steps
CALL rk4(xnh,vnh,tshalf,xn,vn,t+tshalf,w,g,q)
LET d1=abs(xn-xnew)
LET d2=abs(vn-vnew)
LET delta=max(d1,d2)
IF delta<eps then
IF t>tmin then
LET tnew=t+tstep
LET s=s+vnew*tstep
IF Poincare=1 then
LET w1=mod(phi-w*t,2*pi)
LET w2=mod(w*tnew-phi,2*pi)
IF w1<w*tstep then
IF w2<w*tstep then
LET ts=w1/w
CALL rk4(x,v,ts,xp,vp,t,w,g,q)
IF abs(xp)>pi then LET xp=xp-2*pi*abs(xp)/xp
CALL GRAPHPOINT(xp,vp,1)
IF sv=1 then WRITE #2:xp,vp
END IF
END IF
END IF
END IF
LET x=xnew
LET v=vnew
LET t=t+tstep !Expand step size
LET tstep=tstep*.95*(eps/delta)^.25
IF abs(x)>pi then !bring theta back in range
LET x=x-2*pi*abs(x)/x
END IF
ELSE !else reduce step size
LET tstep=tstep*.95*(eps/delta)^.2
END IF
IF t>tmax then
LET average=s/(tmax-tmin)
IF average>0 then
CALL GRAPHPOINT(XINT,VINT,4)
IF sv=1 then WRITE #1:xint,vint
END IF
LET i=1000001
END IF
NEXT i
NEXT vint
NEXT xint
CALL addlegend("g="&str$(g)&" q="&str$(q)&" circle=positive",0,1,"white")
CALL drawlegend
GET KEY VARIABLE
CLEAR
print"press <esc> key to finish"
END
!
!
SUB rk4(x,v,tstep,xnew,vnew,t,w,g,q)
DECLARE DEF accel
LET xk1=tstep*v
LET vk1=tstep*accel(x,v,t,w,g,q)
LET xk2=tstep*(v+vk1/2)
LET vk2=tstep*accel(x+xk1/2,v+vk1/2,t+tstep/2,w,g,q)
LET xk3=tstep*(v+vk2/2)
LET vk3=tstep*accel(x+xk2/2,v+vk2/2,t+tstep/2,w,g,q)
LET xk4=tstep*(v+vk3)
LET vk4=tstep*accel(x+xk3,v+vk3,t+tstep,w,g,q)
LET vnew=v+(vk1+2*vk2+2*vk3+vk4)/6
LET xnew=x+(xk1+2*xk2+2*xk3+xk4)/6
END SUB
!
!
DEF accel(x,v,t,w,g,q)
LET accel=-sin(x)-(1/q)*v+g*cos(w*t)
END DEF
!
SUB PARAMS(W,EPS,TSTEP)
LET W=.6666666666
LET EPS=1E-6
LET TSTEP=.5
END SUB